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Abstract 

The spin-1 Heisenberg model on a triangular lattice with the ferromagnetic nearest, 
J\ = —(1 — p)J, J > 0, and antiferromagnetic third-nearest-neighbor, J3 = pj, 
exchange interactions is studied in the range of the parameter ^ p ^ 1. Mori's 
projection operator technique is used as a method, which retains the rotation sym- 
metry of spin components and does not anticipate any magnetic ordering. For zero 
temperature several phase transitions are observed. At p ~ 0.2 the ground state is 
transformed from the ferromagnetic spin structure into a disordered state, which in 
its turn is changed to an antiferromagnetic long-range ordered state with the incom- 
mensurate ordering vector Q = Q' ~ (1.16,0) at p ~ 0.31. With the further growth 
of p the ordering vector moves along the line Q' — Q c to the commensurate point 
Q c = (^,0), which is reached at p = 1. The final state with an antiferromagnetic 
long-range order can be conceived as four interpenetrating sublattices with the 120° 
spin structure on each of them. Obtained results are used for interpretation of the 
incommensurate magnetic ordering observed in NiGa2S4. 
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This work was motivated by the recent synthesis of crystals NiGa2S4 [T] and Ba3NiSb20g 
[2|3] . These compounds demonstrate interesting magnetic properties, which are mainly de- 
termined by a two-dimensional triangular lattice of Ni 2+ ions with spin S = 1. Both systems 
are characterized by a spin disorder at low temperature and strong antiferromagnetic interac- 
tions (the Curie- Weiss temperature is negative). In particular, in the compound NiGa 2 S 4 the 
experiment on neutron scattering revealed the incommensurate short-range order [1] with 
the in-plane correlation length equal to 6.9 lattice spacings at T = 1.5K. The grain size of 
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the sample was several orders of magnitude larger than the correlation length. To describe 
the observed order the J1-J3 classic Heisenberg model with the dominating third- nearest- 
neighbor (TNN) antiferromagnetic (J 3 ) and weak nearest-neighbor (NN) ferromagnetic in- 
teraction (Ji) was proposed pQ. The ratio J1/J3 can be fitted such that the incommensurate 
ordering vector Q c i of the model coincides with the observed momentum Q eX p of the peak 
in the scattering intensity. The vector Q c i is close to the commensurate ordering vector in 
the case when the NN interaction vanishes. To estimate the values of the exchange constants 
in NiGa2S4 ab initio density functional calculations were performed pTJ. It was shown that 
the second-nearest-neighbor exchange constant is negligibly small, while the TNN exchange 
constant J 3 is anomalously large. The NN exchange constant is smaller than J 3 and possi- 
bly ferromagnetic. Hence the ratio of parameters of the classical J1-J3 model proposed in 
Ref. [1] agrees with numerical estimations of Ref. [I]. However, due to quantum effects the 
classical model cannot give a comprehensive description for the system of S = 1 spins. It is 
known that the magnetic phase diagram of the quantum Heisenberg model with competitive 
interactions is much richer. In particular, the phase diagram contains phases with the short 
range order (SRO), which separates the long range ordered (LRO) phases [5.6.7J (the clas- 
sical phase diagram has only LRO states). Earlier the quantum J1-J3 model was considered 
in the Schwinger-Boson mean field approach [8], where a fitting to experimental results [lj 
was attempted. 

In this article we consider the S = 1 J1-J3 model on a triangular lattice with the ferro- 
magnetic NN (Ji = —(1 — p) J < 0, J > 0) and the antiferromagnetic TNN (J 3 = pj > 0) 
couplings. In the following we use J as the unit of energy. The frustration parameter p 
changes from to 1. We use Mori's projection operator technique [9], which retains the 
rotation symmetry of spin components and does not anticipate any magnetic ordering. The 
used method does not require any approximate representations of the spin operators for 
calculating the spin Green's function. In this approach, this function is represented by a 
continued fraction. The elements of the fraction are calculated in a recursive procedure, 
which is similar to Lanczos' orthogonalization [10]. It was found that the magnetic phase 
diagram of the model consists of the LRO ferromagnet at \Ji\ ^> J 3 , a disordered state, 
incommensurate LRO antiferromagnetic phases with ordering vectors located on the line 
Q' — Qc [(Q' ~ (1-16, 0), Q c = (27r/3, 0)] and a combination of four 120° spin structures on 
non-interacting enlarged triangular sublattices. That spin texture corresponds to the case 
when the coupling between nearest neighbors vanishes (Ji = 0,p = 1). Thus in the range 
< p < 1 the system undergoes three phase transitions. We show that the considered model 
can describe the incommensurate short-range order observed in NiGa2S4 at the value of the 
frustration parameter p = 0.82 and at low but finite temperature. 

The Hamiltonian of the model reads 



where s„ and s„ are the components of the spin-1 operators s n , n and m label sites of the 
triangular lattice, a = ±1. The spin-1 operators can be written as s n = J2a=±i ^l 11 ; cr )( n ) 
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and s„ = -\/2(| n , 0)(n, — er| + |n, cr)(n, 0|), where |n, ±1) and |n, 0) are site states with different 
spin projections. As mentioned above, we take into account the NN and TNN interactions, 
Jnm = Ji J2a ^n,m+a + ^3 J2a ^n,m+A with the vectors a and A=2a connecting the NN and 
TNN sites. Hereafter we use the lattice spacing a = |a| as the unit of length. 

The retarded Green's function reads 

D(kt) = -ie(t)([st(t),s z _ k }), (2) 

where s£ = N^ 1 ^ 2 J2 n e _ * kn s^, N is the number of sites, = e %m s^e~ %m and the angular 
brackets denote the statistical averaging. 

We exploit Mori's projection operator technique pfTU] for calculating the Fourier transform 
of Kubo's relaxation function, 

oo oo 

KKk))« = / tfe^KKJt, Kia = 0(f) J dt'([sl(t'),s\]). 

— oo t 

The Fourier transform of Green's function (j2j) can be obtained from this relaxation function 
using the relation 

D(\au)=u((sl\s z _ k )) u -(sl,s*_ k ), (3) 

where (A,B) = i J °° dt([A(t), B\). In this approach, ({s^ls^))^ is represented as the contin- 
ued fraction 

KKk)).- = > (4) 
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where the elements E n and V n of the fraction are determined from the recursive procedure 
[A n , H] = E n A n + A n+l + K-i^n-i, E n = ([A n , H\,A{) (A n , A\)~ l , 

(5) 

K-i = {A n ,Al){A n ^Ai_ x )-\ V- 1 = 0, A = s z k , n = 0,1,2,... 

The operators Ai obtained in the course of these calculations form an orthogonal set, 
(A h A]) (x S tj . 

Using procedure (j5J) we get 

E = («ai,a* k )(4,a-k) _1 = K,^J>(4,^k) -1 = 0, A 1 = i£, 
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V = 6J[- (1 -p)d[ 7 (k) - 1] + P C 2a [ 1 (2k) - 1] 
Ei = if -is-k)(**k, -is-*)" 1 = 0, 



where 7(k) = | cos^) + | cos cos f ^-^J in the orthogonal system of coordinates, 

Ci = ( s n ls n+a) an d da = ( S n ls n+A) are the spin correlations on the NN and TNN sites, 
respectively. At this point we interrupt the continued fraction and calculate (s£, s -k)- m thus 
taken approximation V\ oc (A 2 , A\) = 0. From this equation we find 



(\i 2 s z k ,-ist k }) = 36J 2 [- (1 -p)d[ 7 (k) - 1] + pda[ 7 (2k) - 1] 



(6) 



The quantity i 2 sl. in the left-hand side of this equation is a sum of terms of the type sf s„ 1 s 
Following Refs. [TTfT2] . we use the decoupling 
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for the case 1 7^ m, n. Here a is the vertex correction. In contrast to the case 5 = | [13J, 
the terms with 1 = n or 1 = m do not cancel each other completely. For 5 = 1 the residual 
terms read 

Pi = ^= E J im (|1> +1) (1> 0|s m - |m, +1) (m, 0|sf - Sl + |m, 0) (m, +1| + s+ |1, 0) (1, +1|) . 

We neglect these terms in the following calculations taking into account that (Pi) = 0. 

Using this approximation for i 2 's^., from Eq. ([6]) we find (s£, si k ) and from Eqs. ([3]) and (TjO) 

we get 



D(kuj) 
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where 



W 2 =36J 2 a|(l-p) 2 [ 7 (k)-l] 
+ p 2 [ 7 (2k)-l] 



Cl r, 2(1 - a) 
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p(l - p) [1 - 7 (k)](C" - 7 (2k)d) + [1 - 7 (2k)](C" - 7 (k)C 2a ; 



d = J (f + 2(4^^) + 2C X + C 2a ), C" = I (2(4^) + 2(4^^) + d + (4 X ^3a) 
and C 2 = | (| + (4 ls n+4a) + 2da + 2(s+ 1 s n + 2d )) , d = a x + a 2 , r = 2a x + a 2 , where a x = 
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(1,0) and a 2 = f|, 2 ) are ^ ne basis vectors of the triangular lattice. As follows from 
Eq. ([7]), the quantity is the frequency of the spin excitations. From Eq. ([8]) we see that 
this frequency tends to zero when k — > 0. 

To find the parameters a, Cx, C 2 , C' 2 , C 2a and C" in Eqs. (GO) and (jHJ) we use the relation 
connecting the spin correlations with Green's function (J7j) 



v tk( „- m) (p-l)[7(k)-l]Ci+p [7(2^-1]^ / a* . 



Five equations for Cx, C 2 , C 2 , C^a, C", which are derived from Eq. (j9]), and the equation 
(s^ l s m l ) = 4/3, (10) 

which follows from the constraint (s^) = 2, form the closed set for calculating all parameters 
of Eqs. (7) and (8) for a finite temperature. 

First we discuss briefly the ground state of the corresponding classical model [H]. The 
classical ground state is the combination of the spiral configurations 

S n = ucos(Q c i n) + vsin(Q c i n), (11) 



where u and v are the arbitrary orthogonal unit vectors. The ordering vector Q c i is chosen 
from the condition of minimal energy. For small values of the frustration parameter, < 
p < 0.2, the spin system is ferromagnetically ordered with Q c i = (see the inset in Fig. 1, 
in which the length of the ordering vector as a function of p is shown). For larger values of 
the frustration parameter Q c i becomes nonzero. The energy minimum is achieved at six, in 
general case incommensurate, vectors (±fc, 0) , (±|&, ±^/cj and f±§&, T^k) ■ One of these 
vectors coincides with the experimental ordering vector Q exp ~ (1.981,0) at p c i = 0.826. 
At p = 1 the vectors Q c i become commensurate and their length is equal to ^l. Let us 
consider the spin configuration with Q c i = (k,0). The dependencies of the scalar products 
between the classical NN (S • S a ) and TNN (S • S 2a ) spins on p are shown in Fig. 1. 
These quantities are classical analogs of the spin-spin correlators C\ and C 2a - In the pure 
ferromagnetic state (|Q c i| = 0) all spins are codirectional. At p > 0.2 the angle between 
spins So and S 2a grows and spins become opposite in direction at p w 0.45. With further 
decrease of the ferromagnetic coupling J\ and rise of the antiferromagnetic interaction this 
angle diminishes. At p — 1 when J x = four triangular sublattices with lattice spacing 
2a become independent and the ground state is the combination of the 120° structures on 
each sublattice. One can say that the ferromagnetic interaction J x at p ~ 0.45 promotes 
the increase of the antiferromagnetic correlations between S and S 2a and suppresses the 
geometric frustration on the four sublattices. 

Let us now turn to the quantum case. The numerical results were obtained by solving the 
system of equations (9), (10) for the entire range of the frustration parameter < p < 1 and 
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for the temperatures T/J = and 0.2 on lattices up to 216 x 216 sites with periodic boundary 
conditions. For the solution of the mentioned set of six equations we used the Optimization 
toolbox of the Matlab package. If the system has LRO at T = the summation over the 
wave vector in (9) can be divided into the contribution yielded at the ordering vector k = Q, 
which is proportional to the condensation part C, and the fluctuation contribution given by 
other wave vectors [12J. C plays the role of the order parameter. In the case | J\\ ^> J3 the 
ground state is ferromagnetic and the ordering vector Q = (0,0). An additional equation 
for calculating C = Cp is the condition that in the vicinity of the T point does not 
contain terms proportional to k 2 . Solving this set of equations we found that the long range 
ferromagnetic order exists for < p < 0.2. The dependence of Cp on p is shown in Fig. 2. 
In the range < p < 0.2 the ferromagnetic condensation part is practically constant, and 
it vanishes abruptly at p m 0.2. For T = the spin-spin correlation functions G\ and Ci a 
also do not change in the ferromagnetic region (see Figs. 3 and 4) and it is evident that this 
behavior is analogous to the classical case (Fig. 1). C\ and Ci a for T/J = 0.2 are also shown 
in Figs. 3 and 4, and their dependencies are smoother. The evolution of the zero-temperature 
spin-excitation spectrum with p is shown in Fig. 5. The sections of this dispersion along the 
k x axis are shown in Fig. 6 for different values of p. Notice that in Fig. 5(a), which corre- 
sponds to p — 0, the dispersion is parabolic near the T point. Such a spectrum is typical 
for the ferromagnetic LRO. The flat region near the T point in Fig. 5(b) points to radical 
changes in the spectrum, which occur for p > 0.2. 

The frequency of magnetic excitations for T > and p > 0.2 vanishes at the T point, 
and it has a minimum at the k-vector, which is located on the line Q' — Q c (hereafter we 
indicate only one of six symmetric directions in the Brillouin zone). The frequency in the 
minimum is small but finite. Such a minimum at p = 0.52 is shown in the inset of Fig. 6 
by the dashed line. The minimum of the Fourier transform of the exchange coupling, which 
defines the classical ordering vector, is also situated on this line. As mentioned above, this 
classical ordering vector is incommensurate for p / 1. It can be supposed that the ground 
state of the quantum model at a large enough antiferromagnetic coupling J 3 also has the 
long-range antiferromagnetic order with an incommensurate ordering vector. To check this 
assumption the following procedure was performed. Since a finite lattice with a discrete set 
of k points was used in our calculation, at first we determined the value of p, at which 
the frequency minimum falls on some allowed incommensurate wave vector Q at T > 0. 
Assuming that for T = the frequency vanishes at this wave vector, the condensation part 
C = Caf wa s calculated in the same manner as in the earlier works [7|T2"] . The result is 
shown in Fig. 2: Caf is nonzero for p > 0.31. The ordering vector corresponding to the 
boundary value p ~ 0.31 is Q = (1.16, 0). Thus, at T = in the range 0.31 < p < 1 the spin 
lattice has incommensurate LRO. Due to the interplay between the geometrical frustration 
and the exchange interactions the dependence of the condensation part on p has a maximum 
at p « 0.75. With increasing frustration parameter the ordering vector moves along the k x 
direction from Q' to the commensurate point Q c = (^-,0j. This wave vector is reached at 
p = 1, as can be seen in Fig. 6. The spectra in Fig. 5(c) and Fig. 5(d) correspond to the 
incommensurate order at p — 0.52 and the commensurate order at p — 1. The length of 
the incommensurate ordering vector Q at p = 0.82 is equal to 1.978, which is close to the 
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experimentally observed length 1.981 of the vector Q eX p- This frustration parameter p is 
also close to the value p c \ obtained in the classical model. The LRO established at p = 1 
can be conceived as four interpenetrating 120° spin structures on sublattices with twice as 
large lattice spacing. Spin orientations on the different sublattices are independent of one 
another. This is seen in Fig. 3 - in this limit C\, the correlation function between spins on 
different sublattices, vanishes. At the same time C2 a tends to the value of spin-spin correla- 
tion between nearest neighbors in the model with nearest-neighbor exchange on a triangular 
lattice [7j (see Fig. 4). 

In addition to the ferromagnetic and the incommensurate antiferromagnetic phases the phase 
diagram contains a region without any LRO. As seen in Fig. 2, this SRO phase exists in the 
range 0.2 < p < 0.31. The LRO phases of the quantum and classical Heisenberg models are 
similar. The main difference is that in the quantum model these phases are separated by the 
phase with SRO. Thus, at zero temperature with the variation of the frustration parameter 
the system undergoes the phase transitions from the ferromagnetic LRO to the SRO, from 
the SRO to the incommensurate antiferromagnetic LRO, which continuously transforms into 
the commensurate LRO. 

Above we found that the considered model has the state with the incommensurate order- 
ing vector Q = (1.978,0), which is close to that observed in the crystal NiGa 2 S 4 . However, 
the state we obtained has the LRO at zero temperature, while the experimental result cor- 
responds to the SRO state with the correlation length £ equal to 6.9 lattice spacings at 
T = 1.5K pQ. In accord with the Mermin- Wagner theorem [15] the considered model has no 
LRO at T > 0. We estimated £ for T ^ 2.8K and found that it is an order of magnitude 
larger than the experimental one. In the crystal, the correlation length is influenced by the 
interlayer interaction, which is absent in the model. Besides, the difference between the exper- 
imental and calculated lengths may be connected with imperfections of the crystal structure 
and some other frustrating interactions. Nevertheless, obtained results demonstrate that the 
S — 1 J\- J3 Heisenberg model is able to describe key features of the experimental results in 
NiGa2S4. Notice also that quasilinear behavior near the minima of the spin-excitation dis- 
persion explains the quadratic temperature dependence of the specific heat, observed in the 
crystal pp. By our estimations, the range of the quadratic dependence extends to T / J rs 0.4. 
Besides, the shape of the calculated uniform susceptibility qualitatively reproduces the ex- 
perimental data pp. 

In summary, Mori's projection operator technique was used for investigating the excitation 
spectrum and spin correlations of the two-dimensional 5 = 1 J1-J3 Heisenberg model on a 
triangular lattice, which takes into account the ferromagnetic nearest-neighbor [J\ — (— 1 + 
p) J < 0] and the antiferromagnetic third-nearest-neighbor (J 3 = pj > 0) interactions. The 
character of the ground state depends on the frustration parameter p. In the range < p < 
0.2, when the ferromagnetic coupling Ji is larger than J 3 , the system is ferromagnetically 
ordered. At p 0.2 the phase transition into a state with short-range order occurs. The 
next transition takes place at p ~ 0.31; the system changes to the long-range order with 
the incommensurate ordering vectors Q, which are located on the lines connecting the V 
point and corners of the hexagonal Brillouin zone. With the growth of p these vectors move 
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to the corners and reach them at p = 1. This commensurate state can be conceived as 
four interpenetrating sublattices with the 120° spin structure on each of them. With p — > 1 
the spin correlations between the sublattices subside. The phases with the long-range order 
of the quantum Heisenberg model are similar to those in the analogous classical model. 
Additionally the quantum model has the phase with the short-range order, which separates 
the long-range ferromagnetic and incommensurate antiferromagnetic phases. The model is 
able to describe the state with the incommensurate short-range order observed in NiGa2S<4. 
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Figures 




P 

Fig. 1. The dependencies of the scalar products between the nearest (So • S a ) and third nearest 
neighbor (So • S2 a ) classical spins on the frustration parameter p. The inset shows the dependence 
of the length of the ordering vector Q c i on p. 



10 




11 




0.0 0.2 0.4 0.6 0.8 1.0 
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Fig. 3. The dependence of the spin-spin correlation function between the nearest-neighbor spins 
Ci = ( s n ls n+ a ) on the frustration parameter p at T = (circles, dashed line) and T/J = 0.2 
(squares, solid line). 
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Fig. 4. The dependence of the spin-spin correlation function between the third-nearest-neighbor 
spins Cia = ( S n ls n+A) on t ne frustration parameter p at T = (circles, dashed line) and T/J = 0.2 
(squares, solid line). 



13 




14 



Figure 5 
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Figure 5 
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Figure 5 



d) 




Fig. 5. The dispersion of spin excitations for different values of the frustration parameter p at 
T = 0. Half of the Brillouin zone is shown. The thick solid line on the base plane is the border 
of the Brillouin zone for the triangular lattice. Panel (d) shows for p = 1. In that case the 
Hamiltonian splits into four terms, each of which describes the Heisenberg antiferromagnet on a 
triangular lattice with twice as large period. The dashed line on the base plane shows one half of 
the Brillouin zone for this lattice. The ordering vector for p = 1 is Q c = (27r/3, 0), which lies at 
the corner of this small Brillouin zone. 
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Fig. 6. The dispersion of spin excitations Wk for different values of the frustration parameter p at 
T = along the k x axis. The inset shows at zero temperature (solid line) and at T/J = 0.2 
(dashed line) for p = 0.52. 
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